#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

class Function
{
public:
    virtual double operator()(double _x) = 0;
    virtual double diff(double _x){return 0;};/// 如果函数不需要导函数，默认导函数为0.
    
};

double poly(vector<double> X, double x, int m)/// 用于生成(x-x_0)*...*(x-x_m)的值
{
    double s = 1;
    for (int i = 0; i <= m;i++)
    {
	s = s*(x-X[i]);
    }
    return s;
};

double poly_diff(vector<double> X, double x, int m)/// 用于生成(x-x_0)*...*(x-x_m)导数的值
{
    double d = 0;
    double s = poly(X, x, m);
    for (int i = 0; i <= m; i++)
    {
	d = d + s/(x-X[i]);
    }
    return d;
};

/** 
 * 根据输入的向量X和函数f生成对应的向量Y，Y[i] = f(X[i])
 * 
 * @param X 
 * @param f 
 * 
 * @return 
 */
vector<double> getvalues(vector<double> X, Function &f)
{
    int n = X.size() - 1;
    vector<double> Y;
    Y.resize(n+1);
    for (int i = 0; i <= n; i++)
    {
	Y[i] = f(X[i]);
    };
    return Y;
};

/** 
 * 牛顿插值的实现
 * 
 * @param X 
 * @param Y 
 * 
 * @return 
 */
class Newton
{
private:
    vector<double> X,Y;
    vector<vector<double>> f;
public:
    Newton(vector<double>& X, vector<double>& Y)
	: X(X), Y(Y)
    {
	int n = X.size();
	f.resize(n+1);
	for (int i = 0; i < n; i++)
	{
	    f[i].resize(n-i);
	    f[0][i] = Y[i];
	}
	for (int i = 1; i < n; i++)
	{
	    int k = 0;
	    for (int j = 0; j < n-i; j++)
	    {
		if (X[j] == X[j+i])
		{
		    f[i][j] = f[0][j+i];
		    k++;
		}
		else
		{
		    f[i][j] = 1.0*(f[i-1][j+1]-f[i-1][j-k])/(X[j+i]-X[j]);
		    k = 0;
		}
	    }
	}
    };
    double operator ()(const double x) const
    {
	int n = X.size() - 1;
	double p = 0;
	for (int i = -1; i < n; i++)
	{
	    p = p+f[i+1][0]*poly(X, x, i);
	}
	return p;
    };
    double diff(double x)
	{
	    double y = 0;
	    int n = X.size() - 1;
	    for (int i = -1; i < n; i++)
	    {
		y = y+f[i+1][0]*poly_diff(X, x, i);
	    }
	    return y;
	};
    void printF()
	{
	    double f0;
	    cout << "f[x_0], f[x_0, x_1], ..., f[x_0, ..., x_n] = " << endl;
	    for (int i = 0; i < f.size(); i++)
	    {
		for (int j = 0; j < f[i].size(); j++)
		{
		    f0 = f[i][0];
		}
		cout << f0 << " ";
	    }
	    cout << endl;
	    cout << endl;
	};
};

